##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)
library(ggsci)
library(patchwork)

##########################################################################################

option_list <- list(
    make_option(c("--gene"), type = "character") ,
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    gene <- "BMP6"
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.molType.IM.tsv",sep="")
    images_path <- paste(work_dir,"/images/mutRatePlot",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene <- opt$gene
mut_rate_gene_file <- opt$mut_rate_gene_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(234,106,79,alpha=255,maxColorValue=255),
  rgb(203,24,30,alpha=255,maxColorValue=255),
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )

names(col) <- c("IM" , "IGC" , "DGC" , "GC")
class_compare <- c("GS" , "CIN" , "MSI")

###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
dat_mutRateGene2 <- dat_mutRateGene

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol==gene)

###########################################################################################
## 填补
for( idn in unique(dat_mutRateGene2$id)){
  if( length(which(dat_mutRateGene$id %in% idn)) == 0 ){
    tmp <- subset(dat_mutRateGene2 , id == idn)[1,]
    tmp$Hugo_Symbol <- gene
    tmp$MutNum <- 0
    tmp$MutRate <- 0
    dat_mutRateGene <- rbind( dat_mutRateGene , tmp )
  }
}

###########################################################################################

dat_plot <- dat_mutRateGene
dat_plot$Molecular.subtype <- factor( dat_plot$Molecular.subtype , levels = c("GS" , "CIN" , "MSI") , order = T )
dat_plot$value_text <- paste0( round(dat_plot$MutRate , 2) * 100 , "%") 

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################

dat_plot$p.value = ""
dat_plot$p_text = ""
dat_plot$Class <- dat_plot$Molecular.subtype
#dat_plot <- subset( dat_plot , Hugo_Symbol %in% unique(dat_plot[dat_plot$MutNum > 2,"Hugo_Symbol"]))
dat_sampleNum <- unique(dat_mutRateGene2[,c("id" , "SampleNum" , "Type" , "Molecular.subtype")])

result <- c()
for(type in unique(dat_plot$Type)){

    for( i in 1:(length(class_compare)-1) ){
        class1 <- class_compare[i]
        for( j in (i+1):length(class_compare) ){
            class2 <- class_compare[j]
            tmp_1 <- subset( dat_plot , Type == type & Class %in% class1 )
            tmp_2 <- subset( dat_plot , Type == type & Class %in% class2 )

            if(nrow(tmp_1)==0){
                tmp_1 <- tmp_2
                tmp_1$Class <- class1
                tmp_1$SampleNum <- dat_sampleNum[dat_sampleNum$Type==type & dat_sampleNum$Molecular.subtype %in% class1,"SampleNum"]
                tmp_1$MutNum <- 0
                tmp_1$MutRate <- 0
                tmp_1$value_text <- ""
            }

            if(nrow(tmp_2)==0){
                tmp_2 <- tmp_1
                tmp_2$Class <- class2
                tmp_2$SampleNum <-  dat_sampleNum[dat_sampleNum$Type==type & dat_sampleNum$Molecular.subtype %in% class2,"SampleNum"]
                tmp_2$Class <- class2
                tmp_2$MutNum <- 0
                tmp_2$MutRate <- 0
                tmp_2$value_text <- ""
            }

            tmp <- rbind( tmp_1 , tmp_2 )
            tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
            p <- fisher.test(tmp_fisher)$p.value

            if( p < 0.001 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
            }
           
            tmp$p.value <- p
            tmp$p_text <- p_text

            result <- rbind( result , tmp )
        }
    }
}


###########################################################################################

barplot_pval <- function(genename,metadata,out_name){
  ## 这里提取了一下需要的内容 至于选择的顺序不知道这样的处理是否合理

  barplotdata <- metadata[c(1,4,5),c("Class" , "MutRate" , "value_text")]
  pvaldata <- metadata[c(1,4,5),c("Class" , "p.value")]
  
  ##设计横杠的位置
  if (max(barplotdata$MutRate)>0.9){
    maxheight <- 1.1
    UNIT <- 0.1
  }else {
    maxheight <- round(max(barplotdata$MutRate),3)
    spare <- (1-maxheight)
    UNIT <- round(spare/7,3)
  }

  height1 <- maxheight + UNIT
  height2 <- height1 + UNIT
  height3 <- height2 + UNIT
  line <- data.frame(x=c(0.5,0.5,3.5),xend=c(2,3.5,2),
                     y=c(height1,height3,height2),yend=c(height1,height3,height2))
  
  ##p值  字的位置和颜色
  trans <- function(num){
    if(num < 0.001){
        up <- floor(log10(num))
        down <- round(num*10^(-up),2)
        text <- paste("P == ",down," %*% 10","^",up)
    }else{
        text <- paste("P == ",round(num , 2))
    }
    return(text)
  } #最后的效果是在类似expression()中呈现的 "=="表示等于 "%*%"表示✖
  
  TEXTUNIT <- UNIT/1.7
  
  if((pvaldata$p.value[1])<0.05){
    curcol="#CA171E"
  }else{
    curcol="#000000"
  }
  pval12 <- data.frame(x=1.2,y=height1+TEXTUNIT,label=trans(pvaldata$p.value[1]),col=curcol)
  
  if((pvaldata$p.value[2])<0.05){
    curcol="#CA171E"
  }else{
    curcol="#000000"
  }
  pval13 <- data.frame(x=2,y=height3+TEXTUNIT,label=trans(pvaldata$p.value[2]),col=curcol)
  
  if((pvaldata$p.value[3])<0.05){
    curcol="#CA171E"
  }else{
    curcol="#000000"
  }
  pval23 <- data.frame(x=2.8,y=height2+TEXTUNIT,label=trans(pvaldata$p.value[3]),col=curcol)
  
  ##设计竖杠的位置
  UNIT <- UNIT/2
  group12 <- data.frame(x=0.5,xend=0.5,x2=2,xend2=2,y=height1,yend=(height1-UNIT))
  group13 <- data.frame(x=0.5,xend=0.5,x2=3.5,xend2=3.5,y=height3,yend=(height3-UNIT))
  group23 <- data.frame(x=2,xend=2,x2=3.5,xend2=3.5,y=height2,yend=(height2-UNIT))
  
  plot <- ggbarplot(barplotdata, x="Class",y="MutRate",
            label = barplotdata$value_text, fill = "Class" ,
            lab.pos="in") +
    xlab(label='') +
    scale_fill_npg() +
    labs(y = 'Proportions (%)') +
    ggtitle(unique(metadata$Type))  +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=line) +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=group12) +
    geom_segment(aes(x = x2, y = y, xend = xend2, yend = yend),data=group12) + 
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=group13) +
    geom_segment(aes(x = x2, y = y, xend = xend2, yend = yend),data=group13) +
    geom_segment(aes(x = x, y = y, xend = xend, yend = yend),data=group23) +
    geom_segment(aes(x = x2, y = y, xend = xend2, yend = yend),data=group23) +
    geom_text(aes(x=x,y=y,label=label),data=pval12,size=4.2,parse = TRUE,colour = pval12$col) +
    geom_text(aes(x=x,y=y,label=label),data=pval13,size=4.2,parse = TRUE,colour = pval13$col) +
    geom_text(aes(x=x,y=y,label=label),data=pval23,size=4.2,parse = TRUE,colour = pval23$col) +
    theme(
        title =element_text(size=5, face='bold'),
        strip.text.x = element_text(size = 10, colour = "black",face='bold') ,
        axis.text.x = element_text(size = 8,color="black",face='bold') ,
        axis.text.y = element_text(size = 8,color="black",face='bold') ,
        axis.title.y = element_text(size = 10,color="black",face='bold') ,
        #axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1 , size = 8,color="black",face='bold'),
        legend.title = element_blank(),
        legend.text = element_text(size = 8),
        legend.position = "none",
        axis.title =element_text(size = 8),axis.text =element_text(size = 7, color = 'black')
      )

    return(plot)
}

metadata <- result[,-c(1:3)]
metadata <- subset(metadata , Type == "All")
metadata$Class <- paste0(  metadata$Class , "(" , metadata$SampleNum , ")")
metadata$Class <- factor( metadata$Class , levels = sort(unique(metadata$Class) ) , order = T )
 
genename <- gene
p1 <- barplot_pval(genename,metadata,out_name)
out_name <- paste0( images_path , "/MutRate_" , gene , ".IM.molType.pdf" )
ggsave( out_name , p1 , width = 3 , height = 4 )

if(1!=1){
    metadata <- result[,-c(1:3)]
    metadata <- subset(metadata , Type == "IM + DGC")
    metadata$Class <- paste0(  metadata$Class , "(" , metadata$SampleNum , ")")
    metadata$Class <- factor( metadata$Class , levels = sort(unique(metadata$Class) ) , order = T )

    genename <- gene
    p2 <- barplot_pval(genename,metadata,out_name)

    metadata <- result[,-c(1:3)]
    metadata <- subset(metadata , Type == "IM + IGC")
    metadata$Class <- paste0(  metadata$Class , "(" , metadata$SampleNum , ")")
    metadata$Class <- factor( metadata$Class , levels = sort(unique(metadata$Class) ) , order = T )

    genename <- gene
    p3 <- barplot_pval(genename,metadata,out_name)

    plot <- p1 + p3 + p2
    out_name <- paste0( images_path , "/MutRate_" , gene , ".IM.molType.pdf" )
    ggsave( out_name , plot , width = 8 , height = 5 )
}